%!Tex Program = xelatex 
\documentclass[a4paper]{article}
%\documentclass[a4paper]{ctexart}
\usepackage{xltxtra}
%\setmainfont[Mapping=tex-text]{AR PL UMing CN:style=Light}
%\setmainfont[Mapping=tex-text]{AR PL UKai CN:style=Book}
%\setmainfont[Mapping=tex-text]{WenQuanYi Zen Hei:style=Regular}
%\setmainfont[Mapping=tex-text]{WenQuanYi Zen Hei Sharp:style=Regular}
%\setmainfont[Mapping=tex-text]{AR PL KaitiM GB:style=Regular}
%\setmainfont[Mapping=tex-text]{AR PL SungtiL GB:style=Regular}

\usepackage{graphicx}
\usepackage{tikz}
\usepackage{pgfplots} 
\usepackage{fontspec}
\usepackage{cite}
\usepackage{geometry} 
\usepackage{amsmath} 
\usepackage{enumitem}
\usepackage{listings}
\usepackage{ctex}
\usepackage{amssymb}

% 接下来是你的其他代码

% 设置中文字体
\setmainfont{AR PL UMing CN}
\geometry{left=20mm,right=20mm,top=20mm,bottom=20mm}

\title{"step-3"教程程序}
\author{江洲}
\begin{document}
\maketitle
\section{ 有限元方法的基本设定}
    这是第一个我们真正使用有限元来计算的例子。我们将解决一个简单的泊松方程，其边界值为零，但右手边非零：\\
    \begin{align}
        -\Delta u &= f \quad \text{in}\  \Omega \\
        u &= 0 \quad \text{on}\  \partial\Omega
    \end{align}
    如果你学过有限元方法的基础知识,就会知道我们要用有限维近似来近似解u。具体来说,\\
    $\bullet$我们首先需要推导出上述方程的弱形式，\\
    $\bullet$可以通过将试函数$\Phi$左乘方程（左乘是为了方便编程，右乘当然也没问题，但是编程会变得复杂），并在域$\Omega$上积分来获得：\\
    $$-\int_{\Omega}\phi \Delta u=\int_{\Omega}\phi f $$\\
    $$-\int_{\Omega}\phi \Delta u=-\int_{\Omega}\phi \nabla^2 u=-\int_{\Omega}\phi \nabla(\nabla u) $$\\
    $$\int_{\Omega}\nabla\phi \nabla u - \int_{\partial \Omega}\phi \mathbf{n}\cdot \nabla u=\int_{\Omega}\phi f $$\\
    试函数$\Phi$必须满足同样的边界条件（用数学术语来说：它需要来自我们所寻求解的集合的切空间），因此在边界上$\Phi = 0$，因此我们要寻找的弱形式为:
    $$(\nabla \phi,\nabla u)=(\phi,f) $$\\
    然后,问题就变成了从适当的空间(这里是Sobolev空间$H^1$)中找出一个函数u,对所有试函数$\Phi$都成立。\\
    当然，在一般情况下，我们无法在计算机上找到这样的函数，而是寻求一个近似值$u_h(x) = \Sigma_j U_j \Phi_j(x)$
    ，其中$U_j$是我们需要确定的未知膨胀系数（这个问题的"自由度"），而$\Phi_j(x)$是我们将使用的有限元形状函数。为了定义这些形状函数，我们需要以下内容：\\
    $\bullet$一个用来定义形状函数的网格。你已经看到如何在step-1和step-2中生成和操作描述网格的对象。\\
    $\bullet$一个描述我们想在参考单元上使用的形状函数的有限元(在deal.II中,它总是单位间$[0,1]$隔，单位正方形$[0,1]^2$单位立方体$[0,1]^3$，取决于你在哪个空间维度计算）。\\
    在step-2中,我们已经使用了类型的对象FE\_Q<2>,它表示通常的拉格朗日单元(Lagrange elements),通过在支持点(support points)上插值定义形状函数。最简单的是FE\_Q<2>,它使用1次多项式。
    在2d中,这些通常被称为双线性,因为它们在参考单元(reference cell)的两个坐标中都是线性的。(在1d中,它们是线性的,在3d中是三线性的;然而,在deal.II文档中,我们一般不做这种区分，而只是简单地将这些函数称为"线性")\\
    $\bullet$一个DoFHandler对象,以有限元对象提供的参考单元描述为基础,枚举网格上的所有自由度。你也已经在step-2中看到了如何做到这一点。\\
    $\bullet$一个映射，用来说明如何从参考单元上的有限元类定义的形状函数中获得实际单元上的形状函数。\\
    默认情况下,除非你明确说明,否则deal.II将使用(双,三)线性映射，所以在大多数情况下，你不必操心这个步骤。\\
    通过这些步骤，我们现在有了一组函数$\Phi_i$，我们可以定义离散问题的弱形式：找到一个函数$u_h$，即找到上面提到的膨胀系数$U_j$，使得：\\
    $(\nabla \phi_i,\nabla u_h)=(\phi_i,f), \quad i = 0,1,\cdots,N-1 $\\
    一个映射，用来说明如何从参考单元上的有限元类定义的形状函数中获得实际单元上的形状函数。
    默认情况下,除非你明确说明,否则deal.II将使用(双,三)线性映射，所以在大多数情况下，你不必操心这个步骤。

    通过这些步骤，我们现在有了一组函数$\phi_i$，我们可以定义离散问题的弱形式：找到一个函数$u_h$，即找到上面提到的膨胀系数$U_j$，使得：\\
    $$(\nabla \phi_i,\nabla u_h)=(\phi_i,f), \quad i = 0,1,\cdots,N-1 $$\\
    请注意,我们在此遵循惯例,即一切从零开始计数,这在C和C++中很常见。如果我们将上式中的
    表示为$u_h(\mathbf{x})=\sum_j U_j\phi_j(\mathbf{x})$，这个方程可以被改写为一个线性系统:\\ 
    \begin{equation}
        \begin{aligned}
        (\nabla \phi_i,\nabla u_h) &= \left(\nabla \phi_i, \nabla\left[\sum_j U_j \phi_j \right] \right) 
        &= \sum_j (\nabla \phi_i, \nabla[U_j\phi_j]) 
        &= \sum_j (\nabla \phi_i, \nabla \phi_j)U_j
        \end{aligned}
    \end{equation}        
    其中,矩阵A和右手边的F定义为:\\
    \begin{equation}   
        \begin{aligned}
        A_{ij} &= (\nabla \phi_i,\nabla \phi_j) \\ 
        F_i &= (\phi_i,f) 
        \end{aligned}
    \end{equation}    

\section{试函数应该左乘还是右乘？}
    在我们继续描述如何计算这些数量之前，请注意，如果我们从右侧而不是从左侧将原方程乘以一个试函数，那么我们将得到一个线性系统，其形式为\\
    $$U^TA=F^T$$\\
    有一个行向量$F^T$。通过转置这个系统，这当然相当于解决了\\
    $$A^TU=F$$\\
    这里和上面一样，因为$A = A^T$\\
    但一般来说，我们不会采用右乘，为了避免任何形式的混淆，经验表明，只要养成左乘方程的习惯（就像数学文献中经常做的那样），就可以避免一类常见的错误，因为矩阵自动正确，在对理论进行实现时不需要转置。\\

\section{组装矩阵和右手边的向量}
    现在我们知道我们需要什么了（即：可以表示矩阵和向量的对象，以及计算$$A_{ij},F_i$$的方法），接下来看看如何实现这些:\\
    $\bullet$$A$的对象是SparseMatrix类型,而$U$和$F$的对象是Vector类型。我们将在下面的程序中看到哪些类是用来解决线性系统的。\\
    $\bullet$我们需要一种方法来进行积分。在有限元方法中，最常见的是使用正交法，也就是说，积分被每个单元上的一组正交点的加权和所取代。也就是说，我们首先把对域$\Omega$的积分分成对所有单元的积分，即\\
    \begin{equation}     
        \begin{aligned}
        A_{ij}&=(\nabla \phi_i,\nabla \phi_j)=\sum_{k\in \mathbb{T}}\int_K\nabla \phi_i \cdot \nabla \phi_j \\     
        F_i&=(\phi_i,f)=\sum_{k\in \mathbb{T}}\int_K\phi_i f 
        \end{aligned}
    \end{equation} 
    然后用正交法对每个单元的贡献进行近似计算：\\
    \begin{equation}     
        \begin{aligned}
        A_{ij}^K &= \int_K\nabla \phi_i \cdot \nabla \phi_j \approx \sum_{q}\nabla \phi_i(\mathbf{x}_q^K) \cdot \nabla \phi_j(\mathbf{x}_q^K) w_q^K \\     
        F_i^K &= \int_K\phi_i f \approx \sum_q \phi_i(\mathbf{x}_q^K) f(\mathbf{x}_q^K) w_q^K 
        \end{aligned}
    \end{equation} 
    首先，我们需要一种方法来描述正交点的位置$x_q^K$和它们的权重$w_q^K$。它们通常以与形状函数相同的方式从参考单元映射出来,即隐含地使用MappingQ1类,或者,准确地说是Mapping类的派生类之一。参考单元上的位置和权重由派生自Quadrature基类的对象来描述。
    通常,人们选择一个正交公式(即一组点和权重),使正交正好等于矩阵中的积分;这可以实现,因为积分中的所有因子都是多项式,由高斯正交公式完成,在QGauss类中实现。\\
    然后我们需要一些东西来帮助我们计算单元$K$上的$\phi_i(\mathbf{x}_q^K)$。这就是FEValues类的作用:它采用一个有限元对象来描述参考单元上的形函数$\phi$,一个正交对象来描述正交点和权重,以及一个映射对象(或隐含地采用MappingQ1类),并提供单元
    中正交点位置上的形状函数的值及其导数，以及积分所需的各种其他信息。\\
    这个介绍的最后一块是要提到,在得到一个线性系统后,要用迭代求解器进行求解,然后进行后处理:我们用DataOut类创建一个输出文件,然后可以用一个常见的可视化程序进行可视化。\\

\section{求解线性系统}
    对于一个有限元程序来说，我们在这里最终得到的线性系统是比较小的：矩阵的大小为$1089\times1089$，这是由于我们使用的网格是$32\times32$，所以网格中有$33^2=1089$个顶点。在后面的许多教程程序中，矩阵大小在几万到几十万之间的情况并不少见。在任何情况下，
    即使是这里的小系统，其矩阵也比在本科生或大多数研究生课程中通常遇到的要大得多，因此问题出现了，我们如何能解决这样的线性系统。\\
    人们通常学习的第一个解决线性系统的方法是高斯消去法。这个方法的问题是，它需要的运算次数与$N^3$成正比，其中$N$是线性系统中的方程或未知数的数量--更具体地说，运算次数是$\dfrac{2}{3}N^3$，相差无几。在$N=1089$的情况下,这意味着我们必须进行大约8.61亿次操作。
    这是一个相当可行的数字,现代的处理器需要不到0.1秒的时间来完成这个任务。但很明显,这是不可能用在大规模计算中的;如果我们在线性系统中有20倍的方程(也就是20倍的未知数),那么就已经需要1000-10000秒或者一个小时的时间了。让线性系统再大十倍,很明显,我们无法在一台计算机上解决它。\\
    相反,我们在这里要做的是采用1952年的一个想法:共轭梯度法(Conjugate Gradient method),或简称为 "CG"。CG是一个"迭代"求解器，因为它形成了一个向量序列，可以收敛到精确解；事实上，在没有舍入误差的情况下，经过$N$次这样的迭代,如果矩阵是对称和正定的,它就可
    以找到精确解。该方法最初是作为另一种精确求解线性系统的方法而开发的,就像高斯消去法一样,但就其本身而言,它没有什么优势,而且在接下来的几十年中基本上被遗忘了。但是,当计算机变得足够强大,可以解决高斯消去法不能很好解决的问题时(20世纪80年代的某个时候),CG被重新发现,
    因为人们意识到它很适合于大型稀疏系统，就像我们从有限元方法中得到的那些。这是因为(1)它计算的向量收敛于精确解，因此我们实际上不必做所有的$N$次迭代来寻找精确解，只要我们对合理的近似值感到满意；(2)它只需要矩阵-向量乘积，这对稀疏矩阵非常有用，因为根据定义，稀疏矩阵
    只有$\mathcal{O}(N)$个条目，因此可以用$\mathcal{O}(N)$次计算完成矩阵-向量乘积，而对密集矩阵做同样操作则需要$N^2$次。因此，我们有希望用最多$\mathcal{O}(N^2)$操作来解决线性系统，而且在许多情况下会少得多。\\
    这些迭代求解器的一个重要组成部分是，我们需要指定要解决线性系统的容差--实质上，是关于我们对愿意接受近似解的误差的一个声明。线性系统$Ax=b$的精确解$x$的近似解$\tilde{x}$的误差定义为$\Vert x - \tilde{x}\Vert$，但这是一个我们无法计算的量，因为我们不知道精确
    解$x$。相反，我们通常把残差，定义为$\Vert b-A\tilde{x} \Vert=\Vert A(x-\tilde{x})\Vert$，作为一个计算措施。然后我们让迭代求解器计算越来越精确的解$\tilde{x}$，直到$\Vert A(x-\tilde{x})\Vert \leq \tau$。一个实际的问题是$\tau$应该取什么值。
    在大多数应用中，设为\\
    $\tau=10^{-6}\Vert b \Vert \\$
    是一个合理的选择。\\

\section{关于实现}
虽然这是你能用有限元方法解决的最简单的方程，但这个程序显示了大多数有限元程序的基本结构，也是几乎所有下面的程序基本上都会遵循的模板。具体来说，这个程序的主类看起来像这样。\\
\begin{lstlisting}[frame=single]  
class Step3
{
  public:
    Step3 ();
    void run ();
 
  private:
    void make_grid ();
    void setup_system ();
    void assemble_system ();
    void solve ();
    void output_results () const;
 
    Triangulation<2>     triangulation;
    FE_Q<2>              fe;
    DoFHandler<2>        dof_handler;
 
    SparsityPattern      sparsity_pattern;
    SparseMatrix<double> system_matrix;
    Vector<double>       solution;
    Vector<double>       system_rhs;
};
\end{lstlisting}
$\bullet$make\_grid()。这是一个可以称为预处理的函数。顾名思义,它设置了存储triangulation对象。在后面的例子中,它还可以处理边界条件、几何形状等。\\
$\bullet$setup\_system()。在这个函数中,所有的数据结构都被设置为解决问题所需要的。特别是,它将初始化DoFHandler对象并正确确定与线性代数有关的各种对象的大小。这个函数通常与上面的预处理函数分开,因为在一个与时间相关的程序中,每当网格被自适应细化时(我们将在step-6
中看到如何做),它可能至少每隔几个时间步就被调用一次。另一方面，在上面的预处理函数中，设置网格本身只在程序开始时执行一次，因此，它被分离成自己的函数。\\
$\bullet$assembl\_system()。这就是计算矩阵和右手边内容的地方，在上面的介绍中已经详细讨论过了。由于对这个线性系统进行处理在概念上与计算其条目有很大不同，我们将其与下面的函数分开。\\
$\bullet$solve()。这就是我们计算线性系统$AU=F$的解$U$函数。在当前的程序中,这是一个简单的任务,因为矩阵是如此简单,但只要问题不再那么微不足道,它就会成为程序大小的一个重要部分(当你对库有了更多的了解时,就会逐渐意识到,可以参阅step-20、step-22或step-31)。\\
$\bullet$output\_results()。最后，当你计算出一个解决方案后，你可能想用它做一些事情。例如，你可能想以可视化的格式输出它，或者你可能想计算你感兴趣的物理量：例如，热交换器中的热通量、机翼的空气摩擦系数、最大桥梁载荷，或者仅仅是某一点上的数值解的值。因此，这个函数
是对你的解决方案进行后处理的地方。\\


\section{关于类型的说明}
deal.II通过命名空间types中的别名定义了一些整数(integral)类型。(在上一句话中，"积分(integral)"一词被用作形容词,与名词integer相对应。它不应该与代表曲线或曲面下的面积或体积的名词"integral" 相混淆。形容词integral在C++世界中被广泛使用，如 "积分类型"(integral
 type)、"积分常数"(integral constant)等。)特别是,在这个程序中,你会在几个地方看到types::global\_dof\_index:一个整数类型,用来表示自由度的全局索引,即在定义在triangulation之上的DoFHandler对象中特定自由度的索引(相对于特定单元中特定自由度的索引)。对于当前的
 程序(以及几乎所有的教程程序)，你将有几千个到几百万个全局未知数(对于$Q_1$元素)在2d的每个单元上有4个未知数,在3d有8个未知数)。因此,允许为全局DoF指数存储足够大的数字的数据类型是无符号整数unsigned int,因为它允许存储0到略高于40亿的数字(在大多数系统中,整数是32位的)。
 事实上,这就是types::global\_dof\_index的内容。\\
\begin{lstlisting}[frame=single]
using types::global_dof_index = typedef unsigned int
\end{lstlisting}
在更实际的情况下，这种类型的存在意味着在装配过程中，我们创建一个$4*4$的矩阵(在2d中,使用$Q_1$元素)来表示我们当前所处的单元的贡献，然后我们需要将这个矩阵的元素添加到全局（系统）矩阵的适当元素中。为此，我们需要获取当前单元的局部自由度的全局索引，为此我们将始终使用下面这段代码：\\
\begin{lstlisting}[frame=single]
cell->get_dof_indices (local_dof_indices);
\end{lstlisting}
其中local\_dof\_indices被声明为:\\
\begin{lstlisting}[frame=single]
std::vector<types::global_dof_index> local_dof_indices (fe.n_dofs_per_cell());
\end{lstlisting}
这个变量的名字可能有点名不副实--它代表"在当前单元上局部定义的那些自由度的全局指数"——但持有这种信息的变量在整个库中普遍是这样命名的.
\section{总结}
\input{data.txt}  
\begin{figure}[htbp]
    \centering
    \includegraphics[width=0.25\linewidth]{output.eps}
    \caption{\label{fig:图像}This is the image of step3.}
\end{figure}

\end{document}
